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We investigate nonlinear localized modes at light-mass impurities in a one-dimensional, strongly- 
compressed chain of beads under Hertzian contacts. Focusing on the case of one or two such 
"defects," we analyze the problem's linear limit to identify the system eigenfrequencies and the 
linear defect modes. We then examine the bifurcation of nonlinear defect modes from their linear 
- - - counterparts and study their linear stability in detail. We identify intriguing differences between 

the case of impurities in contact and ones that are not in contact. We find that the former bears 
similarities to the single defect case, whereas the latter features symmetry-breaking bifurcations 
with interesting static and dynamic implications. 

Ch 1 PACS numbers: 05.45.Yv, 43.25.+y, 45.70.-n, 46.40.Cd 

3 ' 

hj ; 

(N ' INTRODUCTION 

(N ; 

i— 1 ' In the present work, we investigate nonlinear localized modes resulting from configuration heterogeneity in granular 
crystals. This entails a confluence of three key research themes: intrinsic localization through nonlinearity, wave 
propagation in granular chains, and localization through extrinsic disorder. 

Intrinsic localized modes (ILMs), otherwise known as discrete breathers, have been a central theme for numerous 
theoretical and experimental studies over the past two decades Eljl 1, 0, |f| . The original theoretical proposal of ILMs 
in prototypical settings such as anharmonic nonlinear lattices [g, [7| and the rigorous proof of their existence under 
fairly general conditions 0] motivated numerous studies of such modes in a diverse host of applications, including 
t— ( . optical waveguides and photorefractive crystals |9| , the denaturation of the DNA double strand llOl , micromechanical 
cantilever arrays [HI], nanomechanical resonators 11], superconducting Josephson junctions [HI], as well as Bosc- 
Einstein condensates [3], and electrical lattices [l5|, among many others. 

One-dimensional (ID) granular crystals, consisting of closely-packed chains of elastically interacting particles, have 
drawn considerable attention during the past few years. This broad interest has arisen from the wealth of available 
material types/sizes and the ability to tune their dynamic response to encompass linear, weakly nonlinear, and strongly 
nonlinear regimes H, 17. [Til. Il9j|. Such flexibility makes them perfect candidates for many engineering applications, 



including shock and energy absorbing layers 20 , [5l|, [51, [53] , actuating devices j3lj , and sound scramblers H, IH . 



o 

f^*) . Because of these possibilities, it is crucial to investigate the effects of defects (imhomogeneities, beads with different 
masses, etc.), allowin g th e observation of interesting physical responses such as fragmentation, anomalous reflections, 
and energy trapping [2(3, El [H, 0, 0, [51, [51, [53, H [II El . 



It is well-known from solid-state physics that localized vibrations in linear lattices can arise from extrinsic disorder 
that breaks the discrete translational invariance of a perfect crystal lattice [H, HI] . Such "disorder" is also well known 



to introduce interesting wave-scattering effects 37[. This phenomenology arises in a wide host of physical applications, 
including superconductors (38j |. electron-phonon interactions [39f . light propagation in dielectric super-lattices with 
embedded defect layers [4(J , defect modes in photonic crystals |4l| , and optical waveguide arrays [H, 0, EJ . 

In the present work, we aim to investigate the confluence of the preceding research themes by examining the 
interplay of "disorder", which induces localized modes, and nonlinearity in granular crystals. We notice in passing 
that the interaction of impurities with solitary waves or a continuous oscillatory signal in non-loaded (i.e., without 
precompression) monomer chains has been investigated numerically (53 , [53] as well as experimentally in the recent 
work of [II]. In these studies, localized oscillations result from the presence of an impurity of lighter mass than the 
remaining chain particles, during the interaction of the impurity with either a solitary wave or a continuous oscillatory 
signal. However, these localized oscillations were all transient, fading away as soon as the wave left the vicinity of the 
impurity. 

Here, by contrast, we examine long-lived localized breathing oscillations, which form robust nonlinear localized 
modes (NLMs) induced by the impurities, in strongly-compressed granular chains. We demonstrate that their fre- 
quency depends not only on experimental parameters such as the precompressive force and the constitution (material 
and size) of the impurity bead but also on the inherent nonlinearity of the system (i.e., the amplitude of the oscilla- 
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tions). We provide a detailed bifurcation and dynamical analysis of a monoatomic chain with a single lower-mass bead 
("impurity") and extend our considerations to monoatomic chains with a pair of lower-mass impurities. We show 
that the wave dynamics in the case of nearest neighbor "impurities" , which contains consecutive lower-mass beads 
without other intervening particles, differs substantially from that in the case of larger separations between impurities. 
We focus on the case of next- nearest- neighbor impurities, revealing its rich bifurcation structure in which strongly 
asymmetric branches of solutions emerge through symmetry-breaking. We monitor the dynamical manifestation of 
this bifurcation (and associated instability) and examine how it is affected by a potential initial asymmetry in the 
impurity masses. 

The remainder of our presentation is structured as follows. We start by discussing the general theoretical setup of 
the homogeneous (no impurity) model. We then perform linear and nonlinear analyses, first for monoatomic chains 
with a single impurity and then for monoatomic chains with a pair of impurities. Finally, we summarize our findings 
and present some possible future directions. 



MONOATOMIC GRANULAR CHAIN 

The interaction between two adjacent elastic spheres is well-known to be described by Hertz's law [46| . The relation 
between the force Fq exerted on two identical spheres and the distance So between their centers results from geometric 
effects and is given by the nonlinear relation 

= A5l 12 , (1) 

where 

3(l-^ 2 )' (2) 

R is the radius of the beads, E is the material's elastic (Young's) modulus, and v is the Poisson ratio of the bead 
material. 

The dynamics of a ID chain composed of beads of a single type (i.e^ monoatomic chain) is thus described by the 
following system of coupled nonlinear ordinary differential equations [161 ] : 

Mui = A[6 Q + - Uif_l % - A[S + m - 2 , (3) 

where Ui is the displacement of the ith bead from its equilibrium position in the initially-compressed chain, i € 
{2, • • • , TV — 1}, and M is the mass of the beads. The bracket [s] + of Eq. ([3]) takes the value s if s > and the value 
if s < 0, which signifies that adjacent beads are not in contact. 

In contrast with Ref. [jjj ]. which considered unloaded chains, we investigate strongly precompressed chains, in which 
Fq takes large values. Considering small amplitude displacements in comparison with the initial ones caused by the 
precompression force, namely 

hzlZ^ « ! ( 4) 

one can Taylor expand the forces in a power series, in which case keeping the displacement terms to fourth order leads 
to the approximate ( "if 2 — Kj, — if 4" ) model form 

Mili = Kiiui-i - 2u t +Ui-i) 

+ K3 ((ui+i - Ui) 2 ~ - Ui) 2 ) 

+ K 4 ((u i+ i - Ui) 3 + (m-x - Ui) 3 ) , (5) 



where 



K % = 6 -A5l'\ K, = - 6 -A^'\ Ki = ^A8f. (6) 



The equations of motion ([5]) are an example of the celebrated Fermi-Pasta-Ulam (FPU) model j47l . l48l |49| ] . If we 



restrict our consideration to very small amplitudes and velocities, we can neglect all of the nonlinear terms from 
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the equations of motion and keep only the harmonic (K 2 ) term. The spectral band of the ensuing linear chain has 
an upper cutoff frequency of uj m — \fAK 2 jm, which corresponds to the lattice vibration in which the neighboring 
particles oscillate out of phase. 

One of the remarkable features of gradually introducing nonlinearity is that its interplay with discreteness leads to 
the emergence of localized modes even in the absence of any inhomogeneity. Such ILM solutions are generic features 
of a large class of Hamiltonian lattices (which includes the FPU model) [J]. ILMs have been studied extensively in 
monoatomic FPU chains (Hoj . 

One of the canonical mechanisms for the generation of these nonlinear localized modes, is the modulational insta- 
bility (MI) of the band edge plane wave. A detailed analysis of this instability (bifurcation) has shown that the MI of 
the upper cutoff mode manifests itself when the following inequality holds (see Sec. 4.3 of [3| and references therein) 

3X2^4 - 4Kf > 0. (7) 

Using Eqs. © one can show that the above mentioned inequality doesn't hold in our setting, which, in turn, indicates 
that small amplitude ILMs bifurcating from the upper band mode don't exist in monoatomic granular crystals. 
(The existence of dark discrete breathers or large amplitude DBs remains an interesting open question for future 
investigations). 



MONOATOMIC GRANULAR CHAIN WITH AN IMPURITY 



Model 



Consider a ID monoatomic chain of beads that contains an impurity at the fcth site. Suppose that the impurity 
(the fcth bead) is made of the same material as the other particles but has a different radius. In particular, we choose 
the homogeneous host chain composed of spherical stainless steel beads of non-magnetic, 316 type (which have elastic 
modulus E — 193 GPa, Poisson ratio v — 0.3, and density p — 8027.17 kg/m 3 (EH) of radius R = 4.76 mm and 
a spherical steel impurity bead with some other radius. We treat the radius of the impurity as a free parameter, 
though in most cases we will use the radius r = 0.8R. We also suppose that the granular chain is compressed by an 
experimentally-accessible static force of Fq = 22.25 N. 

The presence of the impurity bead in the chain gives rise to a mass defect and leads to changes in the force 
constants that destroy the discrete translational symmetry of the crystal. Recalling that the precompressive static 
force Fq induces an initial displacement Sq between neighboring spheres of the same diameter and defining Si to be 
the displacement that it induces between neighboring spheres of different diameter, equations of motion ([3]) at sites 
k — 1, k, and k + 1 become 

Muk^x = Ai[6 + u k -2 - Uk-itl 2 - A 2 [5i + - Uk}+ , 

muk = A 2 [5i + u k -i - «fe]+ 2 - A 2 [Si + u k - u k+1 ] 3 _/ 2 , 
Muk+i = A 2 [Si + u k - u k+ ifl 2 - Ai[5 + u k +i - tifc+ 2 ]+ 2 , 

A 2E i9 r A 2E (M" m 

Al - 3(l-x 2 ) ' '~ 3(1-^) ' (8) 



where m and r are the mass and the radius of the impurity bead while M and R denote the mass and the radius of 
the remaining beads. 
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Harmonic Potential Approximation (Linear Analysis) 



To understand the underlying linear spectrum of the problem, we linearize ([8]) about the equilibrium state in the 
presence of precompression. This yields 

mui = C\(ui-\ — 2ui — Ui+i) , i £ {k — 1, k, k + 1} , 
Milk-i = Ci(u k -2 - Ufc-i) - C 2 {u k ~i - u k ) , 

mil k = C 2 (Uk-i - 2u fe - Ufe+i) , 
Mii fc+1 = C 2 (w/c - Ufc+i) - Ci(uk+i - u fc+2 ) , 



Co 
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Seeking stationary solutions with frequency u>, we substitute 



for all j into Eqs. ((9]) and obtain the following eigenvalue problem: 



where 



/ 1 
1 




\0 



0\ 




... § ... 



1'2 



t'fc 



1 Uiv_i 

1/ V WAT / 



V = 



/ 1 -1 

-1 2 -1 






V o o 



(9) 
(10) 



= CiV 



/ «1 \ 



V V N ) 



\ 




.. 

-1 2 -1 
-1 1 / 



(11) 











1 1 I C*2 C*2 

1 1 + Ci Ci 
C = | 2|f -ga o 

1 -ga i+^a _i 

The eigenvalue problem (|11[) determines the spectrum of the extended phonon excitations and of the localized defect 
mode centered at the impurity site. In general, the presence of impurity beads can create two types of vibrational 
modes: 

1. Resonance modes, when the mass of the impurity bead is larger than the mass of the rest (m > M). 

2. Localized modes, in the opposite case (to < M). 

Each resonance mode has a frequency within the range of frequencies that constitute the phonon band of the homo- 
geneous host crystal and has a vibration amplitude that is larger in the vicinity of the impurity bead. Each localized 
mode, on the other hand, has a frequency fi m p that lies above the band of the normal modes frequencies of the 
homogeneous host crystal and, as shown in Fig. ([1]), has a vibration amplitude that is large at the impurity site but 
decreases very rapidly with increasing distance. Reference [45| used multiple-scale analysis to obtain the analytical 
approximation 
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FIG. 1: (Color online) (a) The normal mode frequencies of the crystal in the presence of the single lighter-mass defect. The 
presence of the impurity bead leads to the manifestation of a localized mode (see inset) with frequency above the band of the 
perfect (monoatomic) host crystal, (b) Numerically-obtained frequency (solid curve) of the localized mode (/imp) as a function 
of the radius ratio r/R compared with the analytical prediction (dashed curve). As r/R — > 1, one sees that fi mp — ► f m ~ 20.67 
kHz (the horizontal dashed line) and the deviation between the two curves becomes larger. 



for fimp- It is important to observe that our setting is very different from the one in 45|. In particular, the 
experimental setup in [45| has no precompression, so a travelling envelope moves over the defect site. As discussed in 
(451 ] . this pulse acts as a "local precompression force" as it travels, which makes the system weakly nonlinear locally 
and results in localized oscillations at the impurity site. By contrast, in the present setting, the chain is strongly 
compressed by a static force, which acts globally in a nonlinear fashion and allows localized modes to be maintained 
indefinitely without further external excitations after they are excited initially (by, e.g., an actuator or the impact 
of a striker particle). The above distinction between local and global forces is the key difference that leads to very 
long-lived localized oscillations around the impurity bead. Such oscillations can last arbitrarily long in principle, but 
in laboratory experiments the presence of dissipative effects will eventually result in the attenuation of these localized 
modes [52| . 

As illustrated in Fig. [TJb), we find excellent agreement between the analytical expression of Eq. (jT2J) and the 
frequency of the localized mode obtained by the eigenvalue system (fTTj) up to radii ratios = 0.6 (see also Fig. 5 of 
[jjj]). In particular, for the material parameters and the precompressive force discussed above, an impurity bead of 
radius r — 0.6R (for which Q « 0.216) yields a localized mode with frequency fi mp « 31.76 kHz, whereas Eq. (fT2"]) 
predicts f a ~ 30 kHz. On the other hand, for r = 0.8R (implying that ss 0.512), the eigenvalue system (fTT| gives 
fimp ~ 23.28 kHz and Eq. (fT2")) predicts f a ~ 20.03 kHz. To provide additional context, we remark that the upper 

cutoff frequency of the precompressed homogeneous host crystal is given by f m = w 20.67 kHz. It is clear 

that the analytical expression in Eq. (|12p is expected to be a good approximation only for m <C M . Otherwise, one has 
to use the numerically-obtained frequency fimp- Additionally, as the radius of the impurity bead becomes smaller, 
the difference between the frequency fimp of the impurity-induced localized mode and the upper cutoff frequency 
becomes larger. Put another way fimp — * fm as r/R — > 1 (as shown in Fig. HI )), while the localized mode becomes 
concomitantly more extended. 



Continuation and Stability Analysis 



In the previous section, we examined the linear response of the granular crystal in the presence of a single defect. 
We now study Eq. ([8]) directly, to examine the nonlinear behavior of the chain. From a physical perspective, this 
will allow us to examine the interplay of nonlinearity and "disorder" . When a perfect nonlinear lattice supports 
ILMs, the presence of impurities can drastically change the properties of such localized modes, resulting in interesting 
phenomena such as the presence of asymmetric impurity modes in nonlinear lattices with a light- mass defect [531 ] . or 
the existence of stable, nonlinear, heavy-mass impurity modes [5J|. As we have mentioned, our nonlinear lattice does 
not support ILMs [see the inequality (|7|)], so the aforementioned phenomena are not expected to be present. However, 
the linear localization at the defect in conjunction with nonlinearity can result in the presence of robust NLMs. 

More specifically, it is important to consider whether the nonlinearity of the chain can support the existence of 
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localized modes with frequencies f m < f < fi mp . In the linear limit, the chain does not support vibrations with 
frequencies in this regime. To answer that question, we perform a parameter continuation starting from the linear 
localized mode and systematically changing (in small steps) the frequency from f imp towards the upper cutoff f m . 
For each of the intermediate frequencies, we identify NLMs to high precision, via a Newton method in phase space, 
using free boundary conditions and chains with N — 79 beads. In order to identify the relevant branch of solutions, 
we use as an initial guess the localized impurity-induced mode (see insert of Fig. [Ha)), as this was obtained from 
the linear eigenvalue problem (fTTj) . The momenta of all the sites can be fixed to zero, following 55j, due to the time 
reversibility of the system. For details of this continuation method, see Ref. Q and references therein. 

We show the results of the continuation in Fig. [21 which allowed us to obtain localized solutions for all frequencies 
/ G [f m , fimp)- In insets of Fig. [2(a), we show three examples of these solutions; they have frequencies f\ = 22.9 
kHz, fi = 21.65 kHz, and f m = 20.67 kHz. We examined the stability of these localized modes by computing their 
Floquet multipliers \j, which describe the behavior of trajectories near the periodic solution. We show the locations 
in the complex plane of the Floquet multipliers for the three NLM profiles in insets of Fig. [2(a). As is well-known, 
if all eigenvalues Xj have unit magnitude, then the localized periodic solution is linearly stable. However, if |Aj| > 1 
for some j, then a perturbation along the corresponding eigenvector ej grows by the factor |Aj| after one complete 
period. In Fig. H(b), we show one period of the spatiotcmporal evolution of the localized mode with frequency f\. In 
Fig. [2(c), we show the absolute value of the maximal eigenvalue, which is associated with the instability growth rate. 

The family of the localized solutions was found to exhibit an oscillatory instability. In general, oscillatory instabilities 
may arise either due to collision of Floquet multipliers associated with two extended eigenvectors, or between ones 
associated with an extended and one localized. In our case, a careful study of the unstable Floquet multipliers 
and the corresponding eigenvectors reveals that the oscillatory instabilities are caused by the collision of extended 
modes belonging to the two arcs of overlapping continuous phonon spectrum of the Floquet matrix [56j . During the 
continuation of the solutions, typically 3 — 5 quadruplets of eigenvalues abandon the unit circle after the collision 
but return to it jointly soon afterwards in parameter space. As discussed in Ref. [56j], the strength of this kind of 
instabilities should depend on the system size (i.e., the number of beads in the chain) and vanishes in the limit of an 
infinite system. The deviations of the unstable eigenvalues from the unit circle are only up to 0.02, and numerical 
integration of the nonlinear impurity modes up to times 100T (where T is their period) reveals their robustness. 

It is relevant to also note that, among all the Floquet multipliers, two pairs are always located at +1 in the complex 
plane. One of them, the so-called phase mode, describes a rotation of the overall phase of the breather, while the 
second one is due to the conservation of the total mechanical momemtum, an additional integral of motion of the 
FPU chains. As one can see in the insets of Fig. [2(a), the spatial profiles of the NLMs have interesting structure. In 
particular, they are characterized by a kink-shaped distortion of the chain, which is caused by the asymmetry in the 
interaction potential (see Sec. 4.1.4 of jij). This asymmetry is evident in the — Ks — K± approximation of the 
model and it arises directly from the fact that K$ ^ 0. The NLM can be thus viewed as a localized vibration which 
induces this kind of distortion into the granular crystal. As one approaches the edge of the phonon band (i.e., as 
fb * fm), the NLMs become more extended and gradually approach their extended (plane wave) linear counterparts 
at the upper band of the linear spectrum. In this limit, the dc distortion and the maximum of the absolute value of 
the associated Floquet mulipliers also increase. 



Excitation of Nonlinear Impurity Modes 

In the previous section we demonstrated that nonlinear localized modes exist in the gap between the band of phonon 
modes of the perfect crystal and the localized impurity mode. This demonstrates that their frequency depends not 
only on the precompressive force and the material parameters of the defect chain [45| (as is the case in the linear limit) , 
but also on the amplitude displacement of the impurity bead — in other words, on the strength of the nonlinearity. 

In this light, we showcase in the present section what we believe is the easiest way to observe these modes. Our 
method is based on the use of a simple localized initial excitation. At time t — 0, we displace the bead at impurity site 
k by an amount that is strong enough to ensure that the nonlinear terms are no longer negligible in the corresponding 
equations of motion. Meanwhile, we keep all of the remaining sites at rest. We then integrate the equations of motion 
([8| using a fourth-order Runge-Kutta numerical scheme. We expect part of the initially-localized energy excitation 
to spread among the other sites. In order to avoid back-scatter of emitted waves, we consider a chain with N = 500 
beads. 

We generated long-lived localized modes using two proof-of-principle simulations. In the first, we considered a 
relatively weak initial displacement of the impurity site, itfc(O) = <5i/10, which we show in Fig. [3(a). In the second, 
we considered a relatively strong initial displacement of the impurity site, Wfe(0) = S±, which we show in Fig. [3(b). In 
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FIG. 2: (Color online) (a) Continuation diagram of the nonlinear impurity modes. Insets: The spatial profiles and the corre- 
sponding locations of the Floquet multipliers A in the complex plane of the localized modes with frequencies /i = 22.9 kHz, 
/2 = 21.65kHz, and f m = 20.67 kHz. (b) One period of the spatiotemporal evolution of the localized mode with frequency 
/i = 22.9 kHz. (c) The maximum of the absolute value of Floquet multipliers as a function of the frequency /;, of the nonlinear 
impurity mode. 

both cases, we used the precompression and material parameters discussed above. Observe that for the strong initial 
displacement, the nonlinearity causes a substantial distortion of the chain. The frequency of the oscillations of the 
impurity site was about 23.25 kHz and max(ufc/5i) ss 0.065 for the first simulation. For the second simulation, we 
observed a frequency of about 22.47 kHz and max(it/c/5i) « 0.5. Both sets of results are in excellent agreement with 
the continuation analysis discussed above, as that gave frequencies of 23.26 kHz and 22.41 kHz for these particular 
values of max(iik / S\) . These simulations therefore clearly illustrate the excitation of the previously analyzed NLMs. 

MONOATOMIC GRANULAR CHAIN WITH TWO IMPURITIES 
Harmonic Potential Approximation (Linear Analysis) 

We now consider a crystal with two impurities, located at sites k and I. As before, we begin with an analysis of the 
linear spectrum. We first examine the case of two identical impurities, for which we consider impurity beads made 
of the same material (stainless steel) as the ones of the host chain (which is again a perfect crystal) but of smaller 
radius (ri = ri = r = 0.8R). Figure [4] shows the frequencies that correspond to localized modes, which we obtained 
using the harmonic approximation, as a function of the distance between the impurities. Interestingly, for this value 
of radius ratio, when the impurity beads are in contact (I — k = 1), the phonon spectrum has just a single localized 
mode with frequency fi m p ~ 25.28 kHz. As shown in inset (a) in the left panel of Fig. [4l the corresponding mode is 
antisymmetric. Studying the phonon spectrum of the case I — k = 1, as a function of the radius ratio r/R, we found 
that for r/R < 0.8 a symmetric mode leaves the phonon band and becomes progressively more localized as the radius 
ratio is decreased. In particular, for r/R — 0.4 the phonon spectrum contains two localized modes, an antisymmetric 
with fi 62.26fc-ffz and a symmetric with /2 ~ 38A8kHz. When / — k > 2, there are two localized modes even for 
r/R = 0.8. In particular, for I — k = 2, the frequency of the symmetric mode [inset (b) in the left panel of Fig. [4] is 
fi 23.98 kHz, and the frequency of the antisymmetric mode [inset (c) in the left panel of Fig. |4] is ji ~ 22.14 kHz. 
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FIG. 3: (Color online) Spatiotemporal evolution of the displacements of the beads with initial conditions (a) Ufc(0) = <5i/10 
and Uj — for all j k and (b) Uk(0) = 5i and Uj = 0. Top insets: The displacements of the impurity bead (which lies 
at the kih site). Bottom inset: magnifications of the spatiotemporal evolution near the impurity site. Observe that after an 
initial (transient) stage of energy shedding in the form of sound waves, a nonlinear localized breathing mode forms in the 
neighborhood of the defect. 
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FIG. 4: (Color online) Left Panel: The frequencies of the localized modes generated by two identical impurities as a function 
of the distance (number of sites) between them. The insets show corresponding localized modes for (a) I — k = 1, for which 
the impurities are in contact and (b,c) I — k = 2. As the impurities are placed increasingly far apart, we see that /i — — > 0. 
Right Panel: The normal mode frequencies of a granular crystal with two impurities, with I — k — 2, of radii n = 0.775i? and 
r 2 = 0.8 Ji. 



As the two impurities are placed farther apart, j\ — j% — * 0, and fi,f%—* f% m p ~ 23.28 kHz, which is the frequency 
of a single impurity localized mode (see the discussion in the previous section). 

We now consider the phonon spectrum for the case of two different impurities with separation I — k = 2 and bead 
radii r\ — 0.775R and T2 = Q.8R. (As before, the impurity beads are made of the same material as those in the host 
chain.) As one can see by comparing the insets in the left and right panels of Fig. [H even a slight difference in the 
radii of impurity beads results in an asymmetric modification of the corresponding localized modes. In this case, we 
also obtain slightly larger mode frequencies: /i w 24.45 kHz and f% w 22.44 kHz. 



Continuation And Stability Analysis For Two Identical Impurities 

We show the results of our parameter continuation for the case of two identical impurity beads in contact (i — k = 1), 
with radius ration r/R = 0.8 in Fig. We find essentially the same phenomenology as we obtained for granular 
crystals with a single impurity. That is, we obtain a family of weakly-oscillatory unstable localized solutions due to 
finite size effects (the magnitudes of the deviations of the unstable eigenvalues from the unit circle are smaller than 
0.025). Again as before, the modes become wider and the characteristic dc distortion of the chain becomes larger 
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FIG. 5: (Color online) Left panel: Continuation diagram of the nonlinear impurity modes for the case of two identical impurity 
beads in contact (I — k = 1). One set of insets shows the profiles of the localized modes with frequencies /i = 24.75 kHz, 
f% = 22.55 kHz, and f m = 20.65 kHz, and the other set shows the corresponding locations of their Floquet multipliers A in the 
complex plane. Right panel: The maximum of the absolute value of Floquet multipliers as a function of the frequency fb of 
the nonlinear impurity mode. 



as one approaches the frequency f m . We show three examples of this family of localized solutions (with frequencies 
fi = 24.75 kHz, $2 = 22.55 kHz, and f m = 20.65 kHz) in insets of Fig. [5ja). In the rest of the insets, we show the 
locations of their corresponding Floquet multipliers in the complex plane. 

Now consider a granular crystal with two impurities that are not in contact. More specifically, we focus on the 
prototypical case of I — k = 2 and r/R = 0.8. As indicated above, the corresponding phonon spectrum contains two 
localized modes: a symmetric one, shown in inset (b) of the left panel of Fig. |4l and an antisymmetric one, shown in 
inset (c) of the left panel of Fig. [4] In Fig. [6l we show continuation diagrams, which display the frequencies of NLMs 
as a function of the maximum displacement of one of the impurity beads (I = 41, in a chain with N = 79 beads) 
normalized to the characteristic (precompression-induced) displacement 8i. First, we examine the family of solutions 
that arises from the symmetric linear mode at /i 23.98 kHz, see Fig. [5] (a). Stability analysis demonstrates the 
presence of a weak oscillatory instability as in the case of a single and two in-contact impurities. Now consider the 
nonlinear localized solutions that bifurcate from the antisymmetric linear mode, for which fi « 22.14 kHz. This 
family of solutions, which corresponds to the branch A\ of the continuation diagram in panel (b) of Fig. [6l is initially 
weakly unstable (due to the finite-size effects discussed previously). At / ~ 21.56kHz, a pair of Floquet multipliers 
leave the phonon bands. The corresponding eigenmodes are symmetric and become progressively localized [57j as 
the frequency decreases. At / w 21.4AkHz, these two localized modes, collide at the (+1,0) point of the unit circle, 
giving rise to a strong instability (called harmonic instability) which is connected to a bifurcation of the corresponding 
NLM. Two new families of solutions (branches Ai and A3) emerge from this bifurcation which, excluding the kink- 
shaped distortion of the system, are symmetric to each other and weakly unstable due to finite-size effects. Thus, this 
bifurcation is somewhat reminiscent of a pitchfork bifurcation [581 ] . In the case of the newly formed branches A 2 and 
A3 past the bifurcation point, and particularly at / » 21.34kHz, the formed localized pair of eigenmodes enters the 
band of eigenvalues associated with extended perturbations giving rise to a new oscillatory instability. This kind of 
instability although size-dependent, in contrast to the oscillatory instability caused by the collision of two extended 
modes, persists even in the limit of an infinite system 59]. 

It is worth noting that the setting of granular chains with two next-nearest-neighbor impurities (i.e., with I — k = 2) 
is reminiscent of double-well configurations in other contexts. For example, both defocusing and focusing nonlinear 
Schrodinger (NLS) equations with double-well potentials are known to exhibit "symmetry breaking" bifurcations like 
the one discussed above [(5(j. (The defocusing case is relevant to the present set ting .) Moreover, these bifurcations 
have even been observed experimentally in both optical [r3l| and atomic systems \&2\ . 

Examining the temporal dynamics of the unstable antisymmetric mode evinces the symmetry-breaking phenomenon. 
To trigger the relevant instability, we use the wave given by the sum of the unstable solution with / « 21.214 kHz and 
the corresponding unstable localized eigenfunction as an initial condition in the full nonlinear equations of motion. 
Its dynamical evolution, which we show in Fig.[7]Ja), reveals the "symmetry breaking" at t s» 0.4 ms. This is followed 
by alternating oscillations between the two impurity sites (A2 and A3 asymmetric modes). As illustrated in Fig.[7jb), 
the dynamic evolution of the weak oscillatory instability of the asymmetric modes is somewhat similar. As pointed 
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FIG. 6: (Color online) (a) Continuation diagram for the nonlinear impurity modes that originate from the symmetric linear 
mode of a granular crystal with two impurities that are separated by one bead (i.e., I — k = 2). One set of insets shows the 
profiles for the localized modes with frequencies / = 23.68 kHz, / = 22.26 kHz, and / = 20.66 kHz. The other set shows the 
locations in the complex plane of their corresponding Floquet multipliers A. (b) Continuation diagram for the antisymmetric 
solution branch (Ai), showing the generation of asymmetric branches (A.2 and .A3) that originate from the bifurcation at 
/ m 21.44 kHz. The insets show the wave profiles (and their Floquet multipliers) of the antisymmetric localized mode before 
the bifurcation (ft = 21.94 kHz), the antisymmetric mode after the bifurcation (/{, = 21.11 kHz), and the asymmetric mode 
with frequency = 21.11 kHz. (c) The maximum of the absolute value of Floquet multipliers as a function of the frequency of 
the nonlinear impurity mode ft for the symmetric branch, (d) Same as (c), but for antisymmetric (Ai) and asymmetric (A2, 
A3) branches. 



out above, such dynamics is reminiscent of theoretical [60| and experimental [611 ] observations of the instability 
manifestation in NLS equations with double well potentials. 



Continuation And Stability Analysis For Two Distinct Impurities 

As we discussed above, we observe a pitchfork-like bifurcation, indicating the emergence of two families of asym- 
metric localized solutions if the two impurities in the granular chain are identical. Employing our analogy with NLS 
equations, this means that the observed dynamics corresponds to that obtained in a symmetric "double well" that 
can be considered as being induced by the identical, next-nearest-neighbor impurities. A natural generalization is 
to consider how the phenomenology we studied above changes when these impurities are distinct. In this case, the 
"double well" becomes asymmetric, and (from bifurcation theory) one expects to see something analogous to what is 
sometimes called an imperfect pitchfork bifurcation [58j]. Such a scenario, which has been observed in NLS equations 
with asymmetric double wells [g3j], involves an asymmetric perturbation of the pitchfork structure, resulting in a 
saddle-node bifurcation and an isolated branch of solutions. 

To observe this breakdown and obtain the associated modified bifurcation picture, we consider the case of two 
distinct impurities with slightly different radii (namely, r\ — 0.775R and T2 — 0.8R) on the prototypical next-nearest- 
neighbor case of I — k = 2. Recall that we showed the normal mode frequencies for this configuration in the right 
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FIG. 7: (Color online) Spatiotemporal evolution of the bead displacements for (a) the harmonically unstable antisymmetric 
localized mode with ft = 21.214 kHz (for which max|A[ sa 1.437) and (b) for the weakly oscillatory unstable asymmetric 
localized mode with the same frequency (for which max|A| » 1.032). 
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FIG. 8: (Color online) (a) Continuation diagram for the nonlinear impurity mode that originates from f% ~ 22.44 kHz. This 
diagram also shows two additional families of solutions that collide at / « 21.1 kHz. The insets show the spatial profiles (and 
associated Floquet multipliers A) of the localized modes at ft = 22.16 kHz and = 20.9 kHz. In the latter case, we show an 
example from each branch, (b) The maximum of the absolute value of Floquet multipliers as a function of the frequency of the 
nonlinear impurity mode ft for the all of the branches. Observe that the isolated branch that stems from the linear impurity 
mode is only weakly unstable in the frequency interval that we probed, whereas the two branches that arise via the saddle-node 
bifurcation and exist for / < 21.1 kHz encompass a strongly unstable family of solutions and a weakly unstable one. 



panel of Fig. 3J We show the continuation and the stability diagrams for the family of solutions originating from 
/2 ~ 22.44 kHz in Fig. [5] (We omit the continuation of the solutions that emerge from the linear mode with frequency 
/i because it resembles the one shown in the left panel of Fig. [51) As suggested above, this amounts to a saddle-node 
bifurcation. The branches of solutions that are analogous to the A\ and A3 branches from Fig. [5] collide at a critical 
value / w 21.1 kHz and disappear. Linear stability analysis shows that one colliding branch consists of strongly 
harmonically unstable solutions and the other consists of weakly oscillatory unstable solutions. The isolated branch, 
which occurs here because we have broken the perfect pitchfork and arises from the frequency /2 of the linear limit, is 
only weakly unstable instead of exhibiting the strong instability we showed in Fig. [S] Once again, this is reminiscent 
of the NLS phenomenology observed in [631 ] . 



CONCLUSIONS AND FUTURE CHALLENGES 



In conclusion, we have investigated the formation of localized modes due to the interplay of nonlinearity and disorder 
(i.e., the presence of defects) in granular crystals. While previous research has been dedicated to the computational 
and experimental identification of such modes 17J, |4fJ , we believe that the present work is the first one to identify such 
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modes in a numerically exact form and offer a systematic analysis of their linear stability. We have argued that these 
localized modes are likely to be absent in monoatomic chains, and we have illustrated that the inclusion of even a single 
defect can produce a linear defect mode whose interplay with nonlinearity generates a full branch of such solutions. 
We have demonstrated that these waves are robust, and it should be straightforward to generate them in experiments 
for a wide range of initial conditions. We extended our single-defect investigation to multi-defect settings, for which 
we examined the prototypical situation of crystals with two defect sites. We identified a much richer phenomenology 
in such situations, as we observed families of both near-symmetric and near-antisymmetric states. We subsequently 
demonstrated that the latter yields additional families of asymmetric solutions through pitchfork-likc bifurcations, 
which we analyzed for defect pairs containing both identical and distinct impurities. 

Although this paper presents some of the first steps towards a systematic understanding and classification of localized 
breathing modes in granular crystals, there are numerous interesting future directions that should be pursued. First, 
one can consider progressively larger numbers of defects. By analogy with NLS settings 64[, we expect this to 
reveal not only a wealth of additional phenomenology (e.g., new families of waves and more complicated bifurcation 
structures) but also implications for how things progress toward the infinite defect limit. If all defects occur only 
as next- nearest- neighbors, we would obtain a perfect diatomic (two species) crystal in this limit. Such a crystal is 
naturally expected to support intrinsic localized mode solutions in the band gap between its acoustic and optical bands. 
Indeed, our preliminary results indicate that such gap breather solutions do indeed arise. The detailed examination 
of such modes, and their higher-dimensional generalizations (including possibly gap vortex modes [65|), are among 
our future goals. 
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